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Realistic modeling of ionic systems necessitates taking explicitly account of many-body effects. 
In molecular dynamics simulations, it is possible to introduce explicitly these effects through the 
use of additional degrees of freedom. Here we present two models: The first one only includes 
dipole polarization effect, while the second also accounts for quadrupole polarization as well as the 
effects of compression and deformation of an ion by its immediate coordination environment. All 
the parameters involved in these models are extracted from first-principles density functional theory 
calculations. This step is routinely done through an extended force-matching procedure, which has 
proven to be very succesfull for molten oxides and molten fluorides. Recent developments based on 
the use of localized orbitals can be used to complement the force-matching procedure by allowing 
for the direct calculations of several parameters such as the individual polarizabilities. 

I. INTRODUCTION 

Ionic liquids are a class of coulombic liquids 1 - which encompasses several types of compounds. They are often 
separated into just two categories, inorganic molten salts and room-temperature ionic liquids, although this may be 
too restrictive. For example, inorganic molten salts contain various families based on monoatomic anions such as 
molten halidcs or molten oxides, but also other ones in which molecular anions are involved such as carbonates and 
nitrates. In the case of room-temperature ionic liquids, organic cations are involved, which opens the way for an 
almost infinite number of potential electrolytes exhibiting a wide range of physical and chemical properties^ 

Ionic liquids play an important role in many fields of chemistry and physics. Apart from the intrinsic interest for 
such systems, they have been studied by computer simulations in many contexts. For example, molten oxides are the 
main components of magmatic melts^— , molten fluorides and chlorides are investigated for their importance in many 
metallurgical processes^— and for their potential use in the nuclear industr y 10 ^ 1 , and room temperature ionic liquids 
for their use in electrochemical storage applications^ 2 - - — or their solvation properties 1 ^—. 

In order to obtain the thermodynamic and transport properties of interest for all these applications, it is necessary 
to simulate large enough samples for a sufficiently long time, which is out of reach of ab initio molecular dynamics 
methods. This is done instead via classical simulations, in which interactions are described by an analytical force field 
derived from a model of the interactions between closed-shell species. Such a models must account not only for the 
classical electrostatic interaction, but also for three interactions arising from the quantum nature of electrons. The 
exchange-repulsion, or van der Waals (VdW) repulsion is a consequence of the Pauli principle, while the dispersion 
(VdW attraction) arises from correlated fluctuations of the electrons. Lastly, the induction term reflects the distortion 
of the electron density in response to electric fields, which are dominated by polarization effects but may also partially 
represent incipient charge transfer associated with bond formation between closed-shell species^ 

In this paper we review some methodological advances which have been proposed in order to include all these 
effects, and particularly the many-body polarization ones, in accurate models of ionic liquids. The critical physical 
factor to appreciate is that the electronic densities of anions, in particular, are strongly affected by their interactions 
with their environment. For example, the oxide ion, which as an isolated ion is unstable by 8 eV with respect to 
autoionization 2 ^, only exists in the condensed phase because of the confining potential exerted on its electrons by its 
neighbour a 21 i 22 . The confining potential compresses and stabilizes the oxide ion and allows us to treat its interactions 
in the condensed phase as those of a closed-shell species, like an inert-gas atom. However a model which is transferable 
from one material to another, or even between different phases of the same material, must recognise that the confining 
potential, and hence the intrinsic properties of the oxide ion will change between them. Furthermore, in a dynamical 
context, fluctuations in the environmental potential due to the thermal motion of the neighbours mean that the ion's 
electron density is also fluctuating and this effect must be included in the interaction model. The effect is less severe 
for halide ions, where the isolated ion is stable, but where the importance of the environmental effect is indicated by 
the observation that the in-crystal polarizability of the fluoride ion in crystalline LiF is 6.2 a.u. compared 16 a.u. for 
the gas-phase io m 21 i 22 One would expect that for the large molecular anions involved in room temperature systems, 
the effect would be weaker still. Our considerations will be directed at the inorganic molten salt systems for which 
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the environmental potential effects arc pronounced, but with the compensating simplification that the complexities of 
molecular shape are avoided. However, we note that similar methods have recently been applied to develop polarizable 
models for room temperature systems.— 

In the first part of the present paper, we provide the expression for the interaction potential in the framework of 
two models with different complexity: the aspherical ion model, necessary for a transferable description of oxides, and 
the simpler polarizable ion model. We then show how all the parameters involved in these models can be determined 
from first-principles density functional theory calculations through an extended force-matching procedure. In a third 
part, we show recent developments based on the use of localized orbitals, which can complement the force-matching 
procedure by allowing for the direct calculations of several parameters such as the individual polarizabilities. 



II. INTERACTION POTENTIALS 



The functional forms introduced below have been proposed previously on the basis of separate examinations of 
each contribution to the ionic interaction energy in a series of well-directed electronic structure calculations on the 
condcnscd-phasei^ir— The potential is best described as the sum of four different components: charge-charge, disper- 
sion, overlap repulsion and polarization. 

^total ^charge ^dispersion ^repulsion _|_ ^/polarization /-^\ 

Two of these terms will attract more our attention, the repulsion and polarization terms. These contributions 
to the overall energy depend on ionic properties like the ionic radius, for the repulsion, or the ionic multipoles, for 
the polarization. As the electronic density of the ions adapt to the environment that act on it through a confining 
potential or local electric fields, these ionic quantities change in a way so as to minimize the overall energy of the 
system. The response to the environment of these additional degrees of freedom such as ionic radius and multipoles 
gives rise to many-body effects in the potential energy of the system. 

The first term corresponds to the electrostatic interaction between two formal charges, 

ychargo = t^L (2) 

i<j 

where q l is the charge of ion i. The use of formal charges underpins the description of the other interactions as 
characteristic of a closed-shell system^! and the expectation that the interaction parameters should be transferable. 
The dispersion component includes dipole-dipole and dipole-quadrupole terms 



^dispersion = _ £ ^( r «)_^_ + /?(r«)^j J (3) 

where C§ (Cg J ) is the dipole-dipole (dipole-quadrupole) dispersion coefficient, and ffi are Tang-Toennies dispersion 
damping functions^ describing the short-range penetration correction to the asymptotic multipole expansion of 
dispersion^ which take the following form: 



.... n lyij r ij\k 
fc=0 

where the b 1 ^ parameter sets the range of the damping effect. The two other terms which include many-body effects 
may differ depending on the systems of interest, as the magnitude of the response to the environment varies from one 
species to another. In the most complex functional form, the so-called "aspherical ion model" (AIM), which has been 
introduced to model oxide materials j 3 - ' 29 ' 30 the overlap repulsion component is given by 



repulsion 
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where 

p» = r v -Sj-SVvi-S™^, (6) 

and summation of repeated indices is implied. Sets O and M are the sets of oxide anions and cations respectively, so 
that the first three summations represent cation-anion, anion-anion and cation-cation short-range repulsions, respec- 
tively. In the cation-anion term, 5a l is a variable which characterizes the deviation of the radius of the oxide anion i 
from its default value, i.e. its "breathing". jV^} is a set of three variables, for each anion, describing the Cartesian 
components of a dipolar distortion of the ion shape. Similarly, {k^} is a set of five independent variables describing 
the corresponding quadrupolar shape distortions. | k | 2 = k%. x + K 2 y + k 2 zz + 2(K xy + K xz + K yz)i an d the two interaction 

tensors Sa = r^/r*^ and S^l = 3r l J / {r l] ) 2 — 8 a p are also used in equation^ The last summations include the 
self-energy terms, that is the energy cost of deforming the charge density of an ion. /3, £ and 7/ arc effective force 
constants determining how difficult it is for a particular ion to be deformed in a spherical, dipolar or quadrupolar 
way. 

The magnitude and orientation of the spherical breathing, dipolar and quadrupolar deformations will be determined 
by minimization of this energy, with respect to these variables: they will therefore depend on the instantaneous 
positions of neighboring ions and consequently change at each timestep in a molecular dynamics run. 

It was assumed throughout that these shape deformations do not affect significantly the anion-anion and cation- 
cation repulsions, which have been represented by simple exponentials as if the ions were effectively spherical in their 
interactions with other ions of the same type. We have also assumed that short range cation distortion effects can be 
neglected in the case of oxide materials. 

The polarization part of the AIM potential includes both dipolar and quadrupolar contributions 



t ^polarization 
V AIM 
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where k\ = k\ = A ^Jiyi C i , = and k\ = j^^tqt- ct l , B l and C l are the dipolc, dipolc-dipole-quadrupole 

and quadrupolc polarizabilitics of ion i. T Q , ( g 7 ... = V Q ,V ) gV 7 • • • (r 1 - 7 ) -1 arc the multipolc interaction tensors^i with the 
superindex indicating the order of the operator. They are computed using the Ewald summation techniquej 31 ' 32 The 
instantaneous values of the dipole and quadrupolc moments are, again, obtained by minimization of this expression. 
The chargc-dipole and chargc-quadrupole cation-anion asymptotic terms are also corrected for penetration effects 3 - 3 - 
at short-range by using Tang-Toennies damping functions 2 ^ 

n (lj-3 ij\k 

9%^) I ,••«- E^T- (10) 

fc=0 

ff3(r") = l-^'£iV" (ID 



with D and Q standing for the dipolar and quadrupolar parts. Note that the bo (6q) parameters, which determine 
the distance at which the overlap of the charge densities begins to affect the induced multipoles, are the same for 
both anions and cations, while the cd (cq) parameters, which are set to unity in the case of the dispersion interaction 
(equation HJ, measure the strength of the ion response to this effect and therefore depend on the identity of the 
ion. The necessity of larger-than-unity values was indicated by ab initio calculations of the induced multipoles in 
distorted crystals.— The short range induction corrections are usually neglected in both anion-anion and cation-cation 
interactions. 

The AIM potential can be seen to contain several (seventeen per oxide anion) additional degrees of freedom (induced 
dipoles and quadrupoles, and ion shape deformations), which describe the state of the electron charge density of the 
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ions. When calculating the forces on the ions in a molecular dynamics simulation, these electronic degrees of freedom 
should have their adiabatic "Born-Oppcnheimer" values, which minimize the total potential energy, for every atomic 
configuration. The values taken by these extra degrees of freedom is then a complicated functional of the atomic 
configurations so that, even if each term in the total energy is written as a sum of individual or pair components, 
the total system energy includes many-body effects. We search for the ground state configuration of the electronic 
degrees of freedom at each time step, using a conjugate gradients routine, i.e. 

where {^} = {6a N , vj? , (12) 

0, where {£ M } = {^>0 (13) 

(14) 

The convergence criteria typically require the energy at successive steps to differ only in the 8th significant figure for 
polarization and in the 10th for repulsion. The dynamics is thus similar to the so-called Born-Oppcnheimer ab initio 
molecular dynamics. Since the ions remained fixed during the conjugate gradient process many of the terms required 
to evaluate the energy do not need to be recalculated at each minimization step. 

In the case of halide systems, we observed that a simpler functional form was sufficient to describe accurately the 
interactions: In that case, we have used the "polarizablc ion model" (PIM), where the repulsion term does not include 
any ion shape deformation term anymore, which corresponds to the simple exponential form: 



^pulsion 
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^pulsion = ^ A ij e ~a'V (15) 

In addition, good representations of the halides are found with the polarization term only including the terms up to 
the dipoles, 
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With these simplifications, the potential now includes only three additional degrees of freedom per ion, associated 
with the induced dipoles, which are calculated in a single minimization procedure: 
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Such interaction potentials have been introduced by our group in the case of molten chloride a 34 ' 35 and fluorides^ - — 
They were also used by Trullas and co-worker to study a series of silver and copper halidesr^— and the PIM model 
now is implemented in the CP2K code.— Going back to the case of oxides, when attention is restricted to a single 
phase or where similar materials are being compared it is often sufficient to neglect the full complexity of the AIM 
model, and we have also successfully used some PIM-type potentials in the case of amorphous GeO a 45 ' 46 and of 
doped zirconia crystals^ - — The potential developed for SiC>2 by Tangney and Scandolo^i on the basis of similar 
force-matching methods to those we describe below is also worth noting. Their model, which performs very nicely in 
reproducing the properties of various phases of Si02 ) 52 ' 53 differs from the PIM one due to the use of partial charges 
for the ions, which hinders its transferability to other silicate compounds. 



III. OBTAINING THE PARAMETERS FROM FORCE-MATCHING AND MULTIPOLE-MATCHING 

In their pioneering work, Tosi and Fumi developed a series of force field parameters for alkali halides in the framework 
of the Born mode h 54 i 55 These parameters were chosen in order to reproduce solid state data. Their potential have then 
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been used in many simulation works involving alkali halidcs in the molten stated - — This type of empirical approach 
has long been the standard procedure for parameterizing interaction potentials. A new source of data for fitting these 
potentials was made possible by the development of ab initio techniques. For example, the most commonly used 
potential for SiC>2 is based on Hartree-Fock calculations of a single H4Si04 cluster and on several crystal properties^ 
Ocffhcr and Elliott used a similar approach to develop an interaction potential for Ge02^ 

The introduction of the Car-Parrinello method^ opened the way to a new class of simulations, namely the ab initio 
molecular dynamics (AIMD)^ In AIMD simulations, may it be of the Car-Parrinello or of the Born-Oppcnheimer 
type, the introduction of an analytical interaction potential is avoided through the use of density functional theory 
(DFT) calculations. This has an appealing "hands free" aspect, but unfortunately these methods cannot yet access 
the time scales necessary to determine many material properties, in particular in ionic liquids transport properties are 
of primary importance. Nevertheless, AIMD simulations are now considered to be very accurate, and they can be used 
as a reference for classical MD simulations. The idea has therefore emerged to develop interaction potentials yielding a 
trajectory that mimics the one which would be obtained by AIMD. Since molecular dynamics is based on an iterative 
integration of Newton's equation of motion, this can be done if the forces generated by the interaction potential are 
the same as the DFT-calculated ones. This approach,— which is now termed force-matching or force-fitting, has been 
used (and extended) for a decade to develop accurate force fields for ionic liquids in the framework of the PIM and of 
the AIM. We will now focus on practical details of these parametcrizations. 

DFT calculations on periodic systems are based on the Kohn-Sham (KS) method to determine the ground state 
electronic wavefunction {<fi } of a given condensed phase configuration containing N ions. This wavefunction is 
obtained by minimization of the Kohn-Sham energy E KS . The force acting on each atom i, with position r 1 , is then 
extracted from: 

^a.DFT - JT- (19) 

a 

The corresponding classical forces for the same condensed phase configuration are easily obtained for a given set of 
parameters from: 

^ytotal 

^a, classical c\ { (^^0 

where the interaction potential could either derive from the PIM or the AIM, and the additional degrees of freedom 
minimize the total energy for this configuration. The idea behind force-fitting consists of finding the set of parameters 
which will minimize the error made in the classical calculation with respect to the DFT one. If this error is minimal, 
the interaction potential can be considered to be of ab initio accuracy. An important concern when dealing with 
ionic systems is that it is often required to perform studies with varying conditions of temperature, pressure and 
compositions. A prerequisite of the interaction potentials therefore is to be transferable from one thermodynamic point 
to another, which can be achieved by determining the interaction potential parameters on a series of representative 
condensed phase configurations instead of one. In practice, this is done by minimizing the following expression: 

1 Nc 1 Nj I 17* T?i I 2 

2 _ 1 \ "* I r DFT r classical I / 01 \ 

XF ~^U w ^h ^wp (21) 

where N c is the number of configurations and Nj is the number of ions in configuration j. Another quantity that can 
also be determined rather straightforwardly from both the ab initio and the classical calculation is the stress tensor 1 . 
It is therefore possible to extend the fitting datasct by minimizing 

i Nc V I TP - TP I 2 

2 _ 1 L^a(i I "oftDFT afi, classical I , 00 x 

Xs-ttz^ nffi i2 y Al ) 

JV c J = 1 I U aj3,DFT I 

In that case relative weights then have to be assigned to \f an< ^ Xs m the minimization procedure. Of course, 
all the parameters of the interaction potential are involved in the calculation of the force and of the stress tensor, so 



In the case of a classical calculation involving additional degrees of freedom, the correct expression for the stress tensor (as well as for the 
heat current) is obtained by exploiting the Hellmann-Feynman theorem and ignoring any explicit derivatives of the additional degrees 
of freedom with respect to the ion positions. 
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that they can all be fitted at this stage. There is nevertheless a risk of interplay between the various terms of the 
potential, so that one term will do the work which should be attributed to the other. This would contribute to a loss 
of transferability of the potential. 

Consequently, we generalize the force-fitting to require that the induced multipoles on our ions in the model 
reproduce the ab initio ones and determine the parameters in V r P° larlzatlon from this requirement. In DFT calculations 
it is possible to calculate the multipole moments of the individual ions thanks to the maximally localized Wannier 
function (MLWF) formalism^ The MLWFs provide a picture of the electron distribution around atoms which is easily 
interpreted from a chemical point of view. They arc determined by unitary transformations of the KS eigenvectors 

N 



m = E U ™ K> ( 23 ) 



m— 1 

where the sum runs over all the occupied KS states {4>n}ne[i,N]i an d the unitary matrix U is determined by iterative 
minimization of the Wannier function spread f2, which is defined as 

1 N 

^ = -7^lE lo gl s "| 2 ; V«=(C|e- < * r »|C) (24) 

^ ' n=l 

when periodic boundary conditions are applied, where a refers to the coordinates axes, x, y, and z. 

A complete theory of electric polarization in crystalline dielectrics has been developed in recent years,— ~— which 
validates the calculation of the dipole moments of single ions or molecules from the center of charge of the subset of 
MLWF which are localized in their vicinity.— ~— The MLWF centers are computed according to^ 

< a = -^(log Sn , a ) (25) 



and the induced dipole moment of a given ion i is defined, in atomic units, as 

>E 



= -2Y,rt (26) 



where r' n is measured with respect to the position of the ion. In ionic systems the attribution of a Wannier center to 
a given ion can be done unambiguously. Several routes have also been proposed to determine the induced quadrupolc 
moments. Aguado et al. have exploited the fact that the MLWFs are well-localised inside the periodic cell and 
evaluate the components of the quadrupole on an ion i from the real-space integral of the charge densities of the 
MLWFs on that ion: 

6 U = - 2 E / dr I CM I 2 (3r«rg - (r'*) 2 <W)/2 (27) 

n£i v cut 

Here the integral runs over the space within a sphere of radius r cut around i, and r n is the distance from r to the 
nucleus position r l calculated with a minimum image convention^ This method is implemented in the CASTEP 
code. Another method consists in determining the second moment of each Wannier function: 

L 2 ( |(C|e-^°e-'^|C)| 2 \ 

(r a rp) n = ^ln |( ^_,.^, — „,,,.„, , ,.„ ( - ) 

The quadrupolc moment of an ion is then given hy^ 



'a/3 ~ ^ 



(3r' a r' - r' 2 5 aP ) n (29) 



where (r' a r'^) n is measured again with respect to the position of the ion. The advantage of this method, which is 
implemented in the CPMD code, is that there is no need to define a cut-off radius around i (apart from the attribution 
of the Wannier center to the ion). Sagui et al. have extended this approach for the calculation of higher multipole 
moments^ 2 - The generalization of the force fitting procedure to multipoles requires the minimization of the functions: 



1 N C 1 N 3 I % i 12 

1 ^ 1 I Mdft A*ci assical I fQf\\ 
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This fitting procedure has now been applied to an important variety of ionic systems. Here we will describe the 
parameterization of the oxide and fluoride series, for which the AIM and PIM frameworks were chosen respectively. 

An important issue which remains is how to include the dispersion terms (Cg 7 and C$ coefficients) in the fitting 
procedure. It is well-known that the dispersion interactions can not be represented accurately by a local or a semilocal 
exchange-correlation functional, such as that provided by the commonly used LDA and GGA functionals^ We have 
therefore decided not to include the dispersion parameters in the fitting procedure, but rather to fix them to some 
given values. The choice of these values is discussed in the next section. Now the second important question is 
whether it is appropriate or not to include these dispersion effects when calculating the classical forces during the 
fitting procedure. Here our choice depends on the nature of the functional used for the DFT calculations. It is well 
known that the LDA functional leads to overbinding effects (although those are not due to dispersion), while GGA 
functionals like PBE 7 ^ or BLY P 75 i 76 lead to underbinding effects. We have therefore chosen to include the dispersion 
effects in the calculation of the classical forces when using the LDA functional only (the effect is included, but the 
corresponding parameters are kept fixed so that they are not fitted during the procedure). 



0.3 




-0 2,1 ' 1 ' 1 ' 1 ' 1 ' ' 

i).2 -0.1 0.1 0.2 0.3 

F (a.u.) 

x, DFT v J 

Figure 1: The quality of the fits of the ^-component of the forces on the ions to those obtained from the ab initio calculations is 
shown for a series of 19 configurations (including high and low pressure polymorphs of AI2O3, MgO, Si02, MgAl2C>4, MgSiO;?, 
CaO, CaAl 2 Si2 08 and two AI2O3 melts). 

In the case of oxide systems, our objective was to develop a potential suitable for simulation of Earth materials 
under crustal to lower mantle conditions of temperature and pressured The most abundant species are those of the 
so-called CMAS system, i.e. the cations Ca 2+ , Mg 2+ , Al 3+ and Si 4+ . In this work, the LDA functional has been 
used for the DFT calculations : A detailed comparison between GGA and LDA based potentials was undertaken for 
AI2O3; this study showed that the description of the melt was much better for the latter p 7 - 7 - The AIM functional form 
was used for the classical interaction potential. We can therefore summarize the fitting procedure this way: 

1. Generation of a series of typical condensed-phase configurations 

2. DFT calculations on each of these configurations: 

(a) Determination of the ground-state wavefunctions, which gives access to the ab initio forces and stress tensor 
components 

(b) Wannier localization, from which the ab initio induced dipoles and quadrupoles components are calculated 

3. Determination of the dispersion coefficients (see next section) 

4. Minimization of \d and XQ with respect to the parameters of the polarization term (^/polarization^ 

5. Minimization of xf and xs with respect to the parameters of the repulsion term (T/ re P u i slon ) 2 



2 Note that all the \i values could be minimized together but the use of a two-step procedure permits to avoid some cancellation of 
errors in the parameterization. 
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Ion pair 
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Ca 2+ -0 2 ~ 


Mg 2+ -0 2 ~ 


Af' + -0 2 - 


Si 4+_ 2- 


A ij 
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a 1J 
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1.5029 
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B ij 
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b ij 




3.5070 


3.9114 
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3.9812 


C ij 




6283.5 


6283.5 


6283.5 
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c ij 




4.2435 


4.2435 


4.2435 


4.2435 


h ij - h ji 




2.0261 


2.2148 


2.2886 


2.1250 


c ij 

C D 




3.9994 


2.8280 


2.3836 


1.5933 


°Q 




1.5297 


1.9300 


2.1318 


1.9566 


r j 
C Q 




1.6301 


1.3317 


1.2508 


1.0592 




44.372 


2.1793 


2.1793 


2.1793 


2.1793 




853.29 


25.305 


25.305 


25.305 


25.305 


°6 — °8 


1.4385 


2.2057 


2.2057 


2.2057 


2.2057 


D 


0.49566 




P 


1.2325 




c 


0.89219 




V 


4.3646 




a 


8.7671 




C 


11.5124 





Table I: Parameters in the repulsive and polarization parts of the potential for the CMAS system. All values are in atomic 
units. The dipole-dipole-quadrupole polarizability was set to 0. 



The quality of the fits of the x-component of the forces on the ions to those obtained from the ab initio calculations 
is shown for a series of 19 configurations of approximately 100 atoms (including high and low pressure polymorphs 
of A1 2 3 , MgO, Si0 2 , MgAl 2 4 , MgSi0 3 , CaO, CaAl 2 Si 2 8 and two A1 2 3 melt) on figure [2 and the parameters 
obtained for the CMAS system with this procedure are summarized in table |TJ 




1 1 i 1 i 1 i 1 i 1 

5 10 15 20 

q/A 1 



Figure 2: Total static structure factors of MgA^CU at T = 2500 K from MD simulation^ 8 , (lines) compared to experimental 
neutron and x-ray diffraction data (circles) at T = 2423 K.— The neutron structure factors are offset by a value of 0.5. 

The potential reproduces well lattice parameters and elastic constants of various oxides and silicates of the CMAS 
system. Typically, the lattice constants from the simulations are up to a few percent lower than the experimental 
values, which is consistent with the use of LDA for the DFT reference calculations. The model also provides reasonable 
predictions for the thermal expansion and the volume compression of major phases of the Earth's mantle in a wide 
pressure range.—. Liquid state densities of oxide and silicate melts are usually within the error bars of the experimental 
data. The model performs somewhat less well for open network-forming structures close to the pure Si0 2 composition. 

Applications of the potential include the investigation of structural and physical properties of solid and liquid oxides 
relevant in geological or technological context. A number of studies were concerned with the stability and mechanisms 
of phase transitions in magnesium silicates (see^S for a summary) . The atomic structure and dynamics of magnesium 
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Ion type 


Polarizability (a 1 ) 


F" 


7.9 


Li+ 





Na+ 


1.0 


K+ 


5.0 


Rb+ 


8.4 


Cs+ 


14.8 


Ca 2 + 


3.1 


Sr 2+ 


5.1 


y3+ 


3.8 


La 3+ 


7.5 


Zr 4+ 


2.9 



Table II: Fitted polarizabilities for the various ions in molten fluorides. All values are in atomic units. 

and calcium aluminatc liquids was studied in combination with neutron and X-ray scattering experiments using 
aerodynamic levitation techniques j£ir— An example of comparison between simulated^ 8 - and experimental^ total 
structure factors of liquid MgAl204 is shown on figure [2] The almost quantitative agreement between measured and 
computed static and dynamic structure factors underlines the high accuracy and transferability of the potential, which 
was optimized using mostly solid state configurations. The structural studies illustrate that molecular simulations 
are essential to interpret the diffraction data of multicomponent liquids since even two different experiments (neutron 
and X-ray diffraction) do not allow the unambiguous determination of some of the anion-cation nearest neighbor 
distances and cation coordination numbers^ Besides the structural model, the simulations provide insight into the 
relation between melt structure and physical properties, such as the viscosity or the self-diffusivities^L^ 8 . The latter 
are important parameters, e.g., to understand the evolution of the early Earth, which may have been covered by 
a magma ocean composed of silicate liquid. Recent simulations of Mg2SiC>4 liquids using the CM AS potential give 
new constraints on the thermodynamic and rheological properties of the melt in the relevant range of pressure and 
temperature. 85 ' 86 A number of ongoing projects make use of the CMAS potential to study solid-solid and solid-liquid 
interfaces, the structure of complex silicate melts in conjunction with X-ray absorption spectroscopy or with the 
interpretation of trace element partioning data between solid and liquid silicates. 



i 1 1 1 1 1 r 




Figure 3: The quality of the fits of the x-component of the forces on the ions to those obtained from the ab initio calculations 
is shown for a series of 25 configurations (pure LiF, NaF, KF, CaF2, SrF2, YF3, LaF3 and ZrF4, L1F-YF3 and LiF-LaF3 
mixtures) . 

In the case of fluoride systems, we are mainly interested in determining physico-chemical properties of melts in- 
cluding many different components. We have therefore built a transferable PIM interaction potential which includes 
Li+, Na+, K+, Rb+, Cs+, Be 2+ , Ca 2 +, Sr 2+ , Y 3 +, La 3 + and Zr 4 + ions. In this work, the PBE functional has been 
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Ion 


pair 


A' 3 


a lj 


b 13 


— U D 


^3 

CD 




F" 


•F" 


282.3 


2.444 




- 


0.0 


0.0 


F- 


-Li+ 


18.8 


1.974 


1 


834 


1.335 




F" 


-Na+ 


44.9 


1.927 


1 


831 


2.5 


0.022 


F" 


-K+ 


138.8 


2.043 


1 


745 


2.5 


-0.31 


F" 


-Rb+ 


151.0 


1.961 


1 


822 


3.5 


-0.44 


F" 


-Cs+ 


151.1 


1.874 


1 


930 


3.4 


0.49 


F" 


-Ca 2 + 


73.0 


1.859 


1 


732 


1.5 


-0.31 


F" 


-Sr 2 + 


82.5 


1.778 


1 


554 


1.331 


-0.33 


F" 




87.4 


1.832 


1 


847 


1.966 


-0.89 


F" 


-La 3+ 


161.6 


1.867 


1 


614 


1.348 


-0.47 


F" 


-Zr 4 + 


62.6 


1.74 


1 


882 


1.886 


-1.0 


M m+ -M m+ 


1.0 


5.0 






0.0 


0.0 


N+ 


-N+ 


5000.0 


3.0 






0.0 


0.0 



Table III: Fitted parameters for the repulsion and polarization terms for molten fluorides, bg = b l g =1.9 for all the pairs. 
M m+ stands for any cation between Li+, Na+, K+, Ca 2+ , Sr 2+ , Y 3+ , La 3+ or Zr 4+ , while N + is Rb+ or Cs+. All values are 
in atomic units. 

used for the DFT calculations (unlike the case of oxide materials, no test has been made using the LDA). We can 
therefore summarize the fitting procedure this way: 

1. Generation of a series of typical condensed-phase configurations 

2. DFT calculations on each of these configurations: 

(a) Determination of the ground-state wavefunctions, which gives access to the ab initio forces components 

(b) Wannier localization, from which the ab initio induced dipolcs components are calculated 

3. Minimization of Xd with respect to the parameters of the polarization term (yp°i arlzatl ° n ) 

4. Minimization of xf with respect to the parameters of the repulsion term (]/ re P ulslon ) 

5. Determination of the dispersion coefficients (see next section) 

Small differences are observed with the oxide materials situation: Due to the smaller number of parameters involved 
in the fitting procedure and the absence of quadrupole polarization effects in the PIM, only xd and xf had to be 
minimized in this case. In addition the dispersion coefficients can be determined after the fitting step because the 
dispersion term is set to when fitting the repulsion (due to the use of a GGA functional in the DFT calculations). 
Since we are mainly interested in systems which are in the liquid state, this procedure was in fact performed twice: 
Firstly, with a series of crystalline configurations, which led to a first set of parameters. And secondly with a scries 
of liquid state configurations of approximately 100 atoms generated from molecular dynamics simulations with the 
parameters of the first step. The quality of the fits of the x-componcnt of the forces on the ions to those obtained 
from the ab initio calculations is shown for a series of 25 configurations (pure LiF, NaF, KF, CaF2, SrF2, YF3, LaF3 
and ZrF<4, LiF-YF3 and LiF-LaF3 mixtures) on figure [31 and the parameters obtained for the fluoride system with 
this procedure are summarized in tables [TT] and IIIII 

A first test of the potential is provided by a simple comparison of the parameters within a series of chemical species. 
For example, in figure 0] we compare the repulsive term of the potential for various M Tn+ -F _ ion pairs. Firstly, in the 
case of the alkaline cations, we observe that the repulsion wall is pushed towards longer distances in a way which is 
coherent with the increase of the ionic radius of the cation in this scries. More interesting is the second test, which 
compares the period 5 elements Rb + , Sr 2+ , Y 3+ and Zr 4+ . All these ions have the same number of electrons and their 
ionic radius therefore do not differ a lot. The repulsion walls are very close one from each other and they become 
shorter-ranged when the cationic charge is increased, which is in agreement with chemical intuition. 

The interest in molten fluorides is mainly due to their involvement in numerous current (aluminum production) 
and proposed (molten salt reactors) technologies. They are often dangerous and corrosive materials, which are liquid 
at high temperatures, so that performing systematic experimental studies is difficult. The development of accurate 
models which are able to predict the physical properties over a wide range of conditions is therefore very valuable. 



11 



2.0 



1.5 



3 



0.5 



a 1o 



. \\\\ 

AW 


■]■■■■ 


■ AW 
AW 




— Li + 




AW 
AW 


- Na + 






— K + 




■ \ \\\\ 


— Rb + 

— Cs + 




\ \ \ \\ 
\ \ \ \\ 


-- F" 






\ \ \ W 




v \ s \SNv 















1.5 2.0 

r(A) 



2.5 




r(A) 



Figure 4: Comparison of the repulsion terms obtained for various cation-fluoride pairs. Top: Alkaline cations Li + , Na + , 
Rb+ and Cs+; bottom: period 5 elements Rb + , Sr 2+ , Y 3+ and Zr 4+ 
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Figure 5: Comparison of the calculated^ and experimental viscosities of LiF-BeF2 mixtures with varying compositions at 
T = 873 K. 



Among all the molten fluorides, an enormous amount of experimental work has been done on LiF-BeF2 mixtures and 
it is by far the best-characterized fluoride melt . 88 i 89 For this reason, it was the ideal testing ground for evaluating 
our interaction potential development approach and we focused first on this system^ We could show that our 
model performed very well in reproducing the polymeric nuoroberyllate speciation and the vibrational properties in 
these melts42 Calculated dynamic properties (electrical conductivity, viscosity) also agree quantitatively with the 
experimental data)2£ To illustrate this, we show a comparison of the calculated 8 -! and experimental^ 8 - viscosities of 
LiF-BcF2 mixtures with varying compositions at T = 873 K on figure [S] The intcrfacial properties were also tested by 
calculating the surface tension from a simulation of the liquid-vapor interfaced In a second step, a series of simulation 
studies was carried on for other molten fluoride systems. This allowed us to predict heat-transport properties of LiF- 
NaF-KF and NaF-ZrF4 mixtures^ 8 , in order to evaluate their suitability to serve as primary coolant in the advanced 
high-temperature reactor concept^ Our simulations were also used to interpret EXAFS and NMR experimental data 
on fluorozirconatc melts 9 - and to calculate electrical conductivities^, diffusion coefficients^ and internal mobilities^. 



IV. BEYOND FORCE-MATCHING: DIRECT CALCULATION OF PARAMETERS 



Although the fitting procedure presented in the previous section has proven to be very successful, in particular in 
the context of the development of transferable potentials, it can sometimes be enhanced through the direct calculation 
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of some of the parameters involved in the models. Among all the parameters, the dipolar polarizability a has an 
evident physical meaning since it measures the relative tendency of the electron cloud of the ion to be distorted from 
its normal shape by an electric field. Although one could think that the polarizability of a molecule or an ion is an 
intrinsic property, which remains unchanged from the gas phase to any condensed phase, this is not the case due to 
the effect of the confining environmental potential in the latte r 21 ' 22 . 

Again, in order to determine the condensed phase polarizabilities from a direct DFT calculation, the MLWF 
formalism provides a very convenient routed When a small external electric field S is applied to the system, the 
linear response may be characterized by an additional field- induced dipolc moment 8fi l on each individual ion. It is 
convenient to think of the applied field as an optical field, in order to distinguish its effect from that of the static 
fields, which are caused by the permanent charge distributions of the molecules, and to think of Sfi 1 as the net induced 
dipole which is oscillating at the optical frequency. For an electronically insulating material, the induced dipole can 
be written in terms of the total (optical frequency) electric field which acts on it, 



(32) 



where the sum runs over all polarizablc ions j =/= i in the system. In this equation, we have introduced the dipole 
polarizability tensor a 1 of ion i, for a particular condensed phase configuration. Also included is the dipole-dipole 
interaction tensor, T^ 2 \ whose components are defined by T^) — V Q V,3^ij; in practice, for a periodic system, it 
will be computed using the Ewald summation technique! 31 ' 32 The first term on the right-hand side of the equation 
represents the direct contribution of the external electric field to the induced dipole; the second-term is the contribution 
of the reradiatcd electric fields due to the dipoles which are induced in all the other ions (j ^ i) in the sample. In 
principle, higher order induced multipoles also contribute to this expansion, but we will ignore them. In a uniform 
external field, the directly induced higher order multipoles on spherical atoms and ions vanish, and even for molecules, 
their effect is expected to be much smaller than that of the dipoles. 

In DFT calculations on periodic systems, the coupling between the external electric field and the electronic system 
is expressed through the macroscopic polarization of the periodically replicated cel l 98 ' 99 and is defined using the Berry 
phase approach of Resta.^i It is then possible to determine the new partial dipole moment for each species in the 
presence of a field, via another localization step. The field induced dipoles arc calculated from the difference between 
the total molecular dipoles in the presence and absence of the field. 

Equation [32] can be inverted to determine the individual electronic polarizabilities for that particular condensed- 
phase configuration. Consider the application of fields £^ a \ along each Cartesian direction a = x,y, z, and denote by 
{Sfi 1 '^} the corresponding values of the induced dipole moments. The total field f l ^ a > at each position r l can be 
obtained from 



^T< 2 > -Sn j '^ (33) 



which is conveniently evaluated from the electric field given by a dipolar Ewald sum. Finally, the polarizability tensor 
of ion i is given by 

a. 1 = (F i y 1 ■ IT (34) 
where F l and II 1 are second-rank three-dimensional tensors defined as 

K = (35) 

Kp = (36) 

The results obtained for a series of molten fluorides and chlorides are summarized in table IIVI Important environ- 
mental effects occur when the nature of the liquid is changed. In the fluoride melts, for example, the polarizability of 
the fluoride anion shifts a lot, passing from 7.8 a.u. in molten LiF (which is very close to the value obtained from the 
force-fitting procedure, this is also the case for the cation polarizabilities) to 11.8 a.u. in CsF. This is due to differences 
in the confining potential, which affects the electron density around a given anion, and originates from both Coulombic 
interactions and the exclusion of electrons from the region occupied by the electron density of the first-neighbor shell 
of cations i 21 ' 22 ' 101 When passing from one cation to another (for example in the series Li + — >Na + ^K + — >Cs + ), two 
effects are then competing: On the one hand, the anion-cation distance increases, which results in a diminution of the 
confining potential, but on the other hand, the volume occupied by the cation electron density also increases, with an 
opposite effect on the confining potential. Here, the observed increase of polarizability with the size of the cation tends 
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to show that the first effect is the most important. Indeed, the value obtained for CsF is approaching the free F~ 
anion polarizability, which is 16 a.u. These effects have also been observed in similar calculations performed on solid 
oxides^ and protic solvent a 102 ' 103 . They have indirectly been used for building quantitative Lewis acidity scales in 
inorganic materials! 104 ' 105 Such a result means that in order to build a completely transferable interaction potential, 
the model should allow the variation of the polarizability in response to the fluctuations of the environments. First 
steps have been made in that direction in the case of MgO onl y 30 ' 106 



System Ion type Polarizability (a 1 ) 



LiF 


F~ 


7.8 




Li+ 


0.3 


NaF 


F" 


9.6 




Na+ 


1.1 


KF 


F" 


10.7 




K+ 


5.5 


CsF 


F" 


11.8 




Cs+ 


16.3 


NaCl 


cr 


25.1 




Na+ 


1.2 


KC1 


cr 


27.0 




K+ 


5.5 



Table IV: Calculated ionic polarizabilities in the liquid phase. All values are in atomic units. 

Since the centers of the MLWFs allow one to derive atomic dipoles in the condensed phase and hence atomic 
polarizabilities, it is tempting to exploit further these localized orbitals, which are often interpreted in terms of non- 
bonding orbitals or lone pairs, to derive the remaining terms of the force field. A systematic procedure has been 
introduced recently for this purpose, in which a series of approximations leads to the attractive and repulsive van der 
Waals interactions from the knowledge of the geometric properties of the MLWEiSi. 

The central idea of this approach consists of assigning a localized electronic density considered as frozen around 
each nucleus, reconstructed from the associated MLWFs. The exponential decay of MLWFs in electronic insulators^ 
suggests modeling them as Slater orbitals determined solely by their spread Sk and center located at r™ determined 
from the localization procedure^. These two properties can be obtained for each type of MLWF from averages over 

^■ 3 I w l 

liquid configurations. Each MLWF contributes to an electronic density pwk = n^ — e~ Kfc l k I with nt = 2 the 

number of electrons in the orbital and «fc = hi order to obtain orientationally averaged interaction potentials, 

one further proceeds to an orientational average of the density around each nucleus i (see figure [6|) , assuming an 
isotropic distribution of the Nw MLWFs centers with the result: 

N w 

p\r) = y UkKk e ~ Kkr > [(1 + K kr> ) sinh(^r<) - * fc r< cosh( Kfc r<)] (37) 
87rr>r< 

where r> = max(|r — r£'| , dk) and r< = min(|r — r™| , dk), with dk = jr^ — r*| the distance of the kth MLWF center 
to the nucleus. The spread and distance to nucleus of Wannier orbitals for molten NaF, NaCl, KF and KC1 are 
summarized in table [Vj 

Figure [5] also compares the electronic density around CI nuclei with the spherically averaged result given by equation 
1371 The agreement with the mean density is very good even at short distances, the scattering of the measured electronic 
densities being mainly due to thermal fluctuations not reproduced by the rigid density approximation. Exponential 
decay of the molecular density is also apparent over more than two decades and valid in the important region where 
the electronic densities of neighbouring ions overlap. For one of the considered species (K + in molten KF and KG), 
the electronic density around the nucleus, derived from the MLWFs, is better described at large distances by a 
nucleus-centered Slater orbital (spread 1.30 a.u. in KF, 1.34 a.u. in KC1) than by equation [37l Equation l37l provides 
a more accurate description near the nucleus, but not further away. 

Once the localized density around each nucleus in condensed phase is known, the repulsive van der Waals interaction, 
which contains a kinetic and an exchange-correlation part, can be derived using two approximations: The above- 
mentioned frozen density ansatz and the use of a kinetic energy functional. The repulsion energy \/ ro P u i slon j s computed 



14 



le-01 




Figure 6: Electronic density around CI ions in molten KC1 from DFT (scattered plot), compared to the analytic ansatz 
(equation I37[) . which involves only the spread and distance to nucleus of the MLWFs. 



System Ion type {£') (tf) 



KF 


F" 


1.35 


0.555 




K+ 


1.36 


0.685 


NaF 


F" 


1.32 


0.560 




Na+ 


0.85 


0.410 


KC1 


cr 


1.92 


0.869 




K+ 


1.36 


0.680 


NaCl 


cr 


1.90 


0.877 




Na+ 


0.86 


0.414 



Table V: Properties of the Wannier orbitals: Average spread S 1 and distance d l to the nucleus. All values are in atomic units. 



at the DFT level from the superposition of the rigid electronic clouds around each ato m 110 i m and the resulting energy 
surface parametrized as a function of atomic coordinates. For each interatomic separation r, we compute 

^■cpulsion^ = T[ps> y^] + J [psjVps] p g d 3 r; (38) 

where p s = p l + fp is the superposition of frozen atomic densities which have their centers separated by a distance 
r, T and e xc are the kinetic energy and exchange correlation functionals, respectively. Any form can be used for the 
kinetic energy functional T and the exchange-correlation energy density e xc , provided that the latter is the one used 
in the starting DFT calculations. However the resulting T/ ro P ulslon depends on their choice. The LLP kinetic energy 
functional^ is generally the most accurat o m i 113 . Note that improvements of the functionals would automatically 
transfer to the resulting force fields. The resulting y re P ulslon generally decays as Ae~ ar in the region corresponding 
to the first neighbour shell, thus justifying the analytical form of the repulsion term in the PIM and the AIM and 
providing the associated parameters. 

As a first example, we applied this strategy to the case of NaF and NaCl melts, starting from DFT calculations 
with the PBE functional and using the LLP kinetic functional to compute the kinetic energy in equation 1381 The 
relevance of this approach can be appreciated on figure [Jj which reports the mean-square relative error on the forces 
Xp defined in the previous section as a function of the Na + -Cl~/F~ repulsion parameters Ae~ ar , all other parameters 
of the force field remaining fixed. The sets obtained by this method stand well within the narrow windows Xf < 1-0; 
thus justifying a posteriori the frozen density approximation for these systems. The corresponding Xf are 0.241 and 
0.363 for NaF and NaCl, respectively. In the latter, less favourable case, the error is not much larger than the value 
obtained with a force field of the same form (PIM) but with all parameters numerically fitted by force-matching 
(Xf = 0.118), and significantly smaller than the one obtained with the popular Tosi-Fumi potential (xp = 1.250)^^. 

Finally, we now turn to the problem of the dispersion interaction. As mentioned previously dispersion forces are not 



15 




a (a.u.) 

Figure 7: Mean-square relative error Xf on the atomic forces compared to DFT (PBE functional) in the NaCl and NaF melts, 
as a function of the Na + -CP/F _ repulsion parameters Ae~ ar . The set obtained by our method (x) stands well within the 
Xf < 1-0 isosurface in both cases. 

easily captured at the DFT level. While the use of non-local functionals seems to be a promising approach;^ their 
computational cost has so far limited their use in AIMD simulations, for which empirical corrections such as the one 
developed by Grimmei^ are preferred. A promising alternative based on a simpler non-local van der Waals density 
functional has recently been proposed *m In the work concerning the oxide and fluoride series we have used Cq and 
Cg J values obtained from a limited set of ab initio (coupled-cluster or M0ller-Plesset) values on in-crystal ionsii& and 
mixing rules to scale values from one material to another. It is also possible to derive systematically approximate 
parameters for the dispersion term defined in equation [3] Silvestrelli showed that MLWF provide again a convenient 
framework to tackle the dispersion effects and introduced the computation a posteriori of the —C§ l /r®[ interaction 
between Wannier centers^, where C| z depends only on the spread of the MLWFs: 

C f = _J_ / dr / dr' vW^W^ (39) 

327T / J\ r \< rc J\r'\<r' c \J PWk M + \f RWl 0') 

where the cut-off radius for each orbital depends on its spread as 

r c = [1.475 -0.866 In S] S. (40) 

These coefficients are obtained from the MLWFs using the expression proposed by Andcrsson et al*^ for the long- 
range interaction between two separated fragments. One can then obtain the dispersion interaction between a pair of 
atoms as the averaged sum over pairs of MLWFs (for k, I from different sites Assuming an isotropic distribution of 
centers around the nuclei i, j at fixed distance dk,i, we obtain to 2 nd leading order: y dls P erslon — _ ^ n=6 g C% /(r lJ ') n 
where the dispersion coefficients are: 

&j = y, c * ( 41 ) 

ci 3 = £ H4 + dl)c$ l (42) 

kei.iej 

We recall that the distances dk,i result from the orbital localization procedure and are not adjustable parameters. This 
approach thus provides two of the parameters entering in equation [31 leaving only the ones in the damping functions 
to be determined independently. 

All the developments introduced in this section can be used in a procedure which can be summarized this way: 

1. Generation of a series of typical condensed-phase configurations 

2. DFT calculations on each of these configurations: 

(a) Determination of the ground-state wavefunctions, which gives access to the ab initio forces components 
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(b) Wannier localization, from which the MLWFs spreads and positions are extracted and the ab initio induced 
dipoles components are calculated 

3. Reconstruction of the electronic density around each ion (following equation I37[) 

4. Calculation of the y r °P ulslon term parameters for each ion pair (following equation I3"51) 

5. DFT calculations on each of these configurations under an applied electric field £ (one calculation for each 
direction), from which new sets of ab initio induced dipoles components are calculated 

6. Calculation of the individual polarizabilities for each ion (following equation I34[) 

7. Minimization of \f and xd with respect to the damping parameters of the polarization term (]/P ol anzation^ on ^ 

8. Determination of the dispersion coefficients for each ion pair (following equations |3"51 and |4"2]) 

This procedure, which leaves fewer parameters to be fitted than a straightforward application of force-fitting, has up 
to now been tested and validated in the simple case of NaCl, NaF, KC1 and KFpiSZ The parameters are summarized 
in table PVTl (the polarizabilities are given in table HV|) . The dynamic (diffusion coefficients, viscosity, thermal and 
electrical conductivity) and static (density) properties of these molten salts were tested against experimental data: 
An excellent agreement was obtained, showing again the predictive character of such first-principles based approaches. 



System Ion pair 


A ij 


a ij 


°6 




b" 




r" 


KF 


F _ -F _ 


144.0 


1.954 


32.7 


100.6 










F - -K+ 


163.3 


2.052 


34.3 


133.1 


2.052 


6.4 


-1.2 




K+-K+ 


109.5 


2.004 


36.0 


168.7 








NaF 


F~-F~ 


162.9 


2.012 


29.4 


92.1 










F~-Na+ 


301.5 


2.456 


5.6 


13.4 


2.456 


9.5 


-0.4 




Na+-Na+ 


599.0 


3.176 


1.3 


2.2 








KC1 


cr-cr 


187.9 


1.601 


333.0 


2516.2 










cr-K+ 


90.2 


1.636 


102.4 


623.6 


1.636 


3.8 


0.4 




K+-K+ 


244.8 


2.183 


35.9 


166.3 








NaCl 


cr-cr 


205.9 


1.626 


298.4 


2294.1 










cr-Na+ 


94.6 


1.770 


14.2 


66.5 


1.770 


3.5 


0.8 




Na+-Na+ 


546.4 


3.1 


1.3 


2.3 









Table VI: Parameters of the PIM interaction potentials obtained for NaCl, NaF, KC1 and KF from the alternative procedure. 
For all the ion pairs, we have taken bg = bg = a 13 . All values are in atomic units. 



V. CONCLUSION 

In conclusion, we have described in this paper two models of different complexity which can be used to describe 
the interactions in ionic liquids. The first one, the polarizable ion model, includes only one many-body effect, namely 
the dipole polarization of the ions. This model was recently implemented in the CP2K code4i The second one, 
the aspherical ion model, not only extends the inclusion of polarization effects up to the quadrupole level, but also 
includes deformations of the anion in a spherical, dipolar or quadrupolar way for the calculation of the short-range 
repulsion. Depending on the system of interest and the quantities that have to be determined, one has to choose the 
most appropriate model. 

These interaction potentials include several parameters which have to be determined for each ion type / ion pair. 
Since the final objective is to predict physico-chemical properties of the materials of interest, no experimental infor- 
mation should be included in the parameterization stage. Here we have summarized the procedure which was used 
for developing PIM potentials for molten fluorides and AIM potentials for molten oxides and consists in performing 
an extended force-fitting using ab initio density functional theory results as reference data. In the last part we show 
how the use of MLWFs can provide alternative ways for calculating some of the parameters. 

Here we have focused on inorganic molten salts. It is worth noticing that an increasing number of studies on 
room-temperature ionic liquids have involved models similar in spirit to the PIM (they mainly differ due to the 
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choice of Lennard- Jones form to represent the van der Waals interactions) ^22rJM. For the moment the force- matching 
procedure was not used in these works, but we have recently extended it to the case of the l-ethyl-3-methylimidazolium 
tetrachloroaluminate ionic liquid.— 
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